/*
 *  neohooke.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 23/05/2011.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "neohooke.h"

NeoHooke::NeoHooke()
:	HyperelasticInvariants(),
	b(0.0),
	typeVol("")
{}

// ***

NeoHooke::NeoHooke(ifstream &in, const string& na, const string & t)
:	HyperelasticInvariants(na,t),
    b(0.0),
	typeVol("")
{
	string keyword;
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("VOL"))
		{
			in >> typeVol; 
			if(!typeVol.compare("mixed")) in >> b;
		}
		else if(!keyword.compare("E"))	in >> E;
		else if(!keyword.compare("NU"))	in >> nu;
        else if(!keyword.compare("G"))	in >> mu;
        else if(!keyword.compare("LAMBDA"))	in >> lambda;
        else if(!keyword.compare("K"))	in >> K;
		else if(!keyword.compare("DENSITY")) in >> density;
	}

	if (E != 0.0 and nu != 0.0)
    {
		mu = E / ( 2.0 * (1.0 + nu) );
        lambda = nu * E / ((1.0 + nu)*(1.0 - 2.0 * nu));
    }
    else if (mu != 0.0 and nu !=0.0 )
    {
    	E = 2.0 * mu * (1.0 + nu);
    	lambda = 2 * mu * nu / (1.0 - 2.0 *  nu);
    }
    else if (mu != 0.0 and lambda !=0.0 )
    {
    	E = mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu);
    	nu = lambda / 2.0 / (lambda + mu);
    }
    else if (E != 0.0 and K !=0.0 )
    {
    	mu = 3.0 * E * K / (9.0 * K - E);
    	lambda = 3.0 * K * ( 3.0 * K - E) / (9.0 * K - E);
    }
    else if (mu != 0.0 and K !=0.0 )
    {
    	lambda = (3.0 * K - 2.0 * mu) / 3.0;
    	E = 9.0 * K * mu / ( 3.0 * K + mu);
    	nu = ( 3.0 * K - 2.0 * mu) / ( 6.0 * K + 2.0 * mu );
    }
    else if (E != 0.0 and mu != 0.0)
	{
		nu = (E - 2.0 * mu) / ( 2.0 * mu );
	    lambda = mu * (E - 2.0 * mu) / (3.0 * mu - E);
	}
	else if (E != 0.0 and lambda !=0.0 )
	{
	    mu = (E - 3.0 * lambda + sqrt(E*E + 9.0 * lambda * lambda + 2.0 * E * lambda)) / 4.0;
	    nu = 2.0 * lambda / (E + lambda + sqrt(E*E + 9.0 * lambda * lambda + 2.0 * E * lambda));
	}
	else if (nu != 0.0 and K !=0.0 )
	{
	    E = 3.0 * K * ( 1.0 - 2.0 * nu);
	    mu = 3.0 * K * ( 1.0 - 2.0 * nu ) / 2.0 / (1.0 + nu);
	    lambda = 3.0 * K * nu / ( 1.0 + nu);
	}
	else if (nu != 0.0 and lambda !=0.0 )
	{
	    mu = lambda * ( 1.0 - 2.0 * nu) / 2.0 / nu;
	    E = lambda * ( 1.0 + nu ) * ( 1.0 - 2.0 * nu ) / nu;
	}
	else if (K != 0.0 and lambda !=0.0 )
	{
	    mu = 2.0 / 3.0 * (K -lambda);
	    E = 9.0 * K * (K -lambda) / ( 3.0 * K - lambda);
	    nu = lambda / (3.0 * K - lambda);
	}

    if (K == 0.0)
    {
        K = lambda + 2.0/3.0 * mu;
    }
    
    cout << "NeoHooke Material constants have been set to:" << endl;
    cout << "\t" << "E= " << E << endl;
    cout << "\t" << "nu= " << nu << endl;
    cout << "\t" << "mu= " << mu << endl;
    cout << "\t" << "lambda= " << lambda << endl;
    cout << "\t" << "K= " << K << endl;

	hypertype ="coupled";
	
	if(!typeVol.compare("cuad") or !typeVol.compare("log") or !typeVol.compare("mixed"))
	{
		hypertype = "decoupled";
	}

	if (!hypertype.compare("decoupled"))
	{
		Svol.setZero(6);
		Siso.setZero(6);
		Ctvol.setZero(6,6);
		Ctiso.setZero(6,6);
	}
}

// *** 

void NeoHooke::strainEnergy()
{
      // Compressible or nearly incompressible neo-Hookean material from Bonet's book (pag. 16).
	if (J > 0.0) 
	{
		if(!hypertype.compare("coupled"))
		{
			W = mu / 2.0 * (I - 3.0) - mu * log(J) + lambda / 2.0 * log(J) * log(J); 
		}
		else
		{
			Wiso = mu / 2.0 * (Im - 3.0);
			
			if(!typeVol.compare("cuad"))	Wvol = K / 2.0 * (J - 1.0) * (J - 1.0);
			if(!typeVol.compare("log"))		Wvol = K / 2.0 * log(J) * log(J);
			if(!typeVol.compare("mixed"))	Wvol = K / (b * b) * (1.0 / pow(J,b) - 1.0 + b *log(J));
			
			W = Wvol + Wiso;
		}
	}
	else
	{
		cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
		Wvol = -0.5;
		Wiso = -0.5;
		W = -1.0;
	}
}

// ***

void NeoHooke::DstrainEnergy()
{
	if (J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W1 = mu / 2.0;
			W2 = 0.0;
			W3 = ( lambda * log(III) - 2.0 * mu) / (4.0 * III);
		}
		else
		{
			Wiso1 = mu / 2.0;
			Wiso2 = 0.0;
			
			if(!typeVol.compare("cuad"))	dWvol = K * (J - 1.0);
			if(!typeVol.compare("log"))		dWvol = K * log(J) / J;
			if(!typeVol.compare("mixed"))	dWvol = K / (b * J) * (1.0 - 1.0 / pow(J,b));
		}
	}
	else
	{
		W1 = -1.0;
		W2 = -1.0;
		W3 = -1.0;
		
		dWvol = -1.0;
		Wiso1 = -1.0;
		Wiso2 = -1.0;
	}
}

// ***

void NeoHooke::DDstrainEnergy()
{
	if (J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W11 = 0.0;
			W12 = 0.0;
			W13 = 0.0;
			W22 = 0.0;
			W23 = 0.0;
			W33 = ( 2.0 * mu + lambda * (1.0 - log(III)) ) / (4.0 * III * III);
		}
		else
		{
			Wiso11 = 0.0;
			Wiso12 = 0.0;
			Wiso22 = 0.0;
			
			if(!typeVol.compare("cuad"))	ddWvol = K; 
			if(!typeVol.compare("log"))		ddWvol = K * (1.0 - log(J)) / (J * J);
			if(!typeVol.compare("mixed"))	ddWvol = K / (b * J * J) * ( (b + 1.0) / pow(J,b) - 1.0);
		}
	}
	else
	{
		W11 = -1.0;
		W12 = -1.0;
		W13 = -1.0;
		W22 = -1.0;
		W23 = -1.0;
		W33 = -1.0;
		
		ddWvol = - 1.0; 
		Wiso11 = - 1.0;
		Wiso12 = - 1.0;
		Wiso22 = - 1.0;
	}
}

// ***

void NeoHooke::DDDstrainEnergy()
{}

//***

//double NeoHooke::getE()
//{	return E;}
//
//// ***
//
//double NeoHooke::getNu()
//{	return nu;}
//
//// ***
//
//double NeoHooke::getLambda()
//{	return lambda;}
//
//// ***
//
//double NeoHooke::getMu()
//{	return mu;}
//
//// ***
//
//double NeoHooke::getDensity()
//{	return density;}
//
//// ***
//
//double NeoHooke::getBulkmodulus()
//{	return K;}
